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Abstract 

In  this  paper  we  develop  a  high  order  explicit  finite  difference  weighted  essentially  non- 
oscillatory  (WENO)  scheme  for  solving  a  hierarchical  size-structured  population  model  with 
nonlinear  growth,  mortality  and  reproduction  rates.  The  main  technical  complication  is 
the  existence  of  global  terms  in  the  coefficient  and  boundary  condition  for  this  model.  We 
carefully  design  approximations  to  these  global  terms  and  boundary  conditions  to  ensure  high 
order  accuracy.  Comparing  with  the  first  order  monotone  and  second  order  total  variation 
bounded  schemes  for  the  same  model,  the  high  order  WENO  scheme  is  more  efficient  and 
can  produce  accurate  results  with  far  fewer  grid  points.  Numerical  examples  including  one  in 
computational  biology  for  the  evolution  of  the  population  of  Gambussia  affinis,  are  presented 
to  illustrate  the  good  performance  of  the  high  order  WENO  scheme. 
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1  Introduction 


In  this  paper  we  develop  a  high  order  finite  difference  weighted  essentially  non-oscillatory 
(WENO)  scheme  for  a  hierarchical  size-structured  population  model  given  by  the  following 
equations 

ut  +  (g(x,  Q(x,  t ))  u)x  +  m(x ,  Q(x,  t ))  u  =  0,  (x,  t )  G  (0,  L]  x  (0,  T] 

g(0,  Q(0,  f))u(0,  f)  =  C(f)  +  f  /3(x,Q(x,t))u(x,t)dx,  f  G  (0,  T]  (1.1) 

io 

u(x,  0)  =  w°(x),  x  G  [0,  L] 

where  u(x,t)  is  the  density  of  individuals  having  size  x  at  time  t,  and  the  non-local  term 
Q(x,t )  is  defined  by 

px  pL 

Q(x,  t)  —  a  u(£)u(£,t)d£+  u(£)u(€,t)d£,  0  <  a  <  1  (1.2) 

Jo  Jx 

for  some  given  function  u.  Q(x,t)  depends  on  the  density  u  in  a  global  way  and  is  usually 
referred  to  as  the  environment.  This  global  dependence  makes  the  design  of  a  high  order 
WENO  scheme  complicated.  Another  complication  is  the  boundary  condition  at  size  x  —  0, 
which  involves  the  function  g  representing  the  growth  rate  of  an  individual,  and  a  global 
dependency  on  the  density  u(x,t )  for  all  x  G  (0,  L\.  The  function  m  in  (1.1)  represents  the 
mortality  rate  of  an  individual.  The  function  (3  in  the  boundary  condition  of  (1.1)  represents 
the  reproduction  rate  of  an  individual,  and  the  function  C  represents  the  inflow  rate  of 
zero-size  individual  from  an  external  source.  We  assume  that  the  functions  g,  m  and  /3  are 
functions  of  both  the  size  x  and  the  environment  Q,  which  in  turn  depends  globally  on  the 
density  u,  hence  the  problem  is  highly  nonlinear. 

The  hierarchical  structured  population  model  (1.1)  describes  population  dynamics  in 
which  the  size  of  an  individual  determines  its  access  to  resources  and  hence  to  its  growth 
or  decay.  This  dependency  is  based  on  the  environment  which  is  a  global  function  of  the 
density  for  all  sizes.  We  refer  to,  e.g.  [5]  for  a  more  detailed  discussion  of  the  background 
and  application  of  the  hierarchical  size-structured  population  models.  These  models  are  used 
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extensively  in  biological  applications.  For  example,  in  [2],  the  model  is  applied  to  study  the 
evolution  of  a  population  of  Gambussia  affinis. 

Hierarchical  structured  population  models  have  been  studied  in  the  literature  in,  e.g. 
[4,  5,  6,  7,  10,  12,  19],  usually  with  more  restrictive  assumptions  on  the  functions  g,  (3  and 
m.  For  example,  in  [5]  the  model  (1.1)  was  considered  for  the  special  situation  g  =  g(Q), 
(3  =  (3(Q),  m  =  m(Q)  and  C  —  0.  In  [4],  the  model  (1.1)  was  studied  with  the  functions  g 
and  (3  depending  linearly  on  the  size  x,  m  independent  of  re,  and  C  =  0.  In  [12],  (1.1)  was 
investigated  with  a  =  0.  The  model  (1.1)  with  the  complete  generality  as  stated  above  was 
studied  in  [1]  and  [14] .  In  [1]  an  implicit  first  order  finite  difference  scheme  was  analyzed  and 
its  stability  and  convergence,  as  well  as  the  existence,  uniqueness  and  well-posedness  (in  L 1 
norm)  of  bounded  variation  weak  solutions  for  (1.1)  were  proved.  In  [14]  we  developed  and 
analyzed  two  explicit  finite  difference  schemes,  namely  a  first  order  upwind  scheme  and  a 
second  order  high  resolution  scheme,  for  solving  (1.1).  Stability  and  convergence  were  proved 
for  both  schemes  in  [14], 

Even  though  the  second  order  high  resolution  scheme  developed  in  [14]  is  significantly 
more  accurate  than  first  order  schemes,  the  performance  can  be  further  improved  by  an 
even  higher  order  scheme.  Of  course,  since  the  solution  of  this  model  can  be  discontinuous, 
a  nonlinearly  stable  high  order  scheme  which  can  maintain  high  order  accuracy  in  smooth 
regions  and  have  sharp,  monotone  discontinuity  transitions  would  be  desirable.  In  this 
paper  we  develop  an  explicit  finite  difference  WENO  scheme  for  solving  (1.1),  following 
the  successful  class  of  WENO  schemes  for  computational  fluid  dynamics  and  for  general 
conservation  laws  [11,  15,  16,  17].  The  main  technical  complication  for  this  development  is 
the  existence  of  global  terms  in  the  coefficient  and  boundary  condition  of  this  model.  We 
carefully  design  approximations  to  these  global  terms  and  boundary  conditions  to  ensure 
high  order  accuracy.  We  then  provide  numerical  examples  to  demonstrate  the  capability 
of  the  scheme  in  solving  smooth  and  discontinuous  solutions.  For  a  smooth  solution,  we 
verify  that  the  complicated  global  boundary  condition  and  coefficients  are  all  implemented 
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accurately  and  the  designed  high  order  accuracy  can  be  achieved.  The  numerical  examples 
also  include  one  in  computational  biology  for  the  evolution  of  the  population  of  Gambussia 
affinis ,  for  which  our  results  with  the  fifth  order  accurate  WENO  scheme  have  a  comparable 
resolution  as  the  second  order  accurate  scheme  used  in  [2]  with  far  fewer  grid  points. 

We  recall  ([1]  and  [14])  that  the  following  assumptions  are  made  on  the  model  functions: 

•  (HI)  g{x,Q)  is  twice  continuously  differentiable  with  respect  to  x  and  Q,  g(x,Q )  >  0 
for  x  G  [0,  L),  g(L,  Q )  =  0,  and  #q(x,  Q)  <  0. 

•  (H2)  m(x,  Q )  is  nonnegative  continuously  differentiable  with  respect  to  x  and  Q. 

•  (H3)  (3{x,Q)  is  nonnegative  continuously  differentiable  with  respect  to  x  and  Q.  Fur¬ 
thermore,  there  is  a  constant  u\  >  0  such  that  sup^Qwg  nxr0)OO\  /3(x,  Q)  <  u>\. 

•  (H4)  uj(x)  is  nonnegative  continuously  differentiable. 

•  (H5)  C(t)  is  nonnegative  continuously  differentiable. 

•  (H6)  u°  E  BV[0,L\  and  u°(x)  >  0. 

In  section  2,  we  present  the  detailed  construction  of  an  explicit,  fifth  order  accurate 
finite  difference  WENO  scheme  for  solving  (1.1).  Section  3  contains  numerical  examples 
demonstrating  the  capability  of  this  WENO  scheme.  Concluding  remarks  are  given  in  section 
4. 

2  A  fifth  order  accurate  finite  difference  WENO  scheme 

First  we  briefly  describe  the  notations  that  we  will  use.  We  assume  the  spatial  domain  [0,  L\ 
is  divided  into  N  cells  with  cell  boundary  points  denoted  by  Xj,  for  0  <  j  <  N,  xo  =  0  and 
Xn  =  L.  For  simplicity  of  presentation  we  will  assume  that  the  mesh  is  uniform  of  size  Ax, 
namely  Xj  =  j  Ax.  Our  scheme  can  also  be  designed  on  a  smoothly  varying  non-uniform 
mesh.  We  denote  the  time  step  by  At.  In  fact,  this  time  step  At  =  A tn  =  tn+l  —  tn  could 
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change  from  one  step  to  the  next,  based  on  stability  conditions,  but  we  use  the  same  notation 
At  without  the  superscript  n  since  we  will  only  consider  one-step  time  discretization  (Runge- 
Kutta).  We  shall  denote  by  u”  and  Q'j  the  finite  difference  approximations  of  u(xj,tn)  and 
Q(xj,tn),  respectively.  We  also  denote 


g?  =  g(xj, Qj)>  m  =  Pin, Q ?),  <  =  m (*;> 9?)»  c"  =  c(f") 


and  we  define  the  standard  discrete  L1  and  L°°  norms  of  the  grid  function  un  by 

N 

Hull  =  |u?|Ax,  ||un||oo  =  max  |<|. 

J  0<j<N  J 

3=1 

For  semi-discrete  approximations  the  superscript  n  referring  to  the  time  level  is  absent. 

We  now  describe  our  development  of  a  fifth  order  finite  difference  WENO  scheme  for 
solving  (1.1),  following  the  construction  of  finite  difference  WENO  schemes  in  [11,  15,  16]. 
Notice  that  we  choose  fifth  order  accuracy  since  it  is  the  most  often  used  WENO  scheme, 
however  our  scheme  can  also  be  designed  for  other  orders  of  accuracy  along  the  same  lines, 
see,  e.g.  [3].  The  semi-discrete  fifth  order  high  finite  difference  WENO  scheme  is  defined  by 
the  following  conservative  approximation 

JtUj  +  zL  (/J+1/2  ~  A-1/2)  +  miui  =  °’  (2-1) 

where  the  fifth  order  accurate  numerical  flux  fj+1/2  is  defined  by  a  weighted  average  of  three 
third  order  accurate  fluxes 


fj+1/2  =  Wifj+i/2  +  w2fj+i/2  +  w3ff+i/2 ,  J  =  °,  1,  •  • '  ,  N- 

The  coefficients  w  1,  W2  and  W3  are  called  the  nonlinear  weights. 

The  three  third  order  accurate  fluxes  fj+1/ 2,  fj+1/2  an(i  fj+1/2  are  obtained  following  the 
procedure  in  [11,  16]  and  are  given  explicitly  by 

f  fj+ 1/2  =  i/i-2  -  If 3- 1  +  T fi 

\  fj+1/2  =  -If 3- 1  +  I  fi  +  If 3+ 1  J  =  2,  •  •  •  ,  IV  -  2, 

{  fj+1/2  =  \f/  +  If 3+1  ~  if 3+  2 
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in  the  middle  of  the  computational  domain.  For  the  two  points  from  the  left  boundary,  these 
fluxes  are  given  by 


f  fl,2  =  \fo  +  lh-\h  [  fl,2  =  ~\fo+  §/l  +  S/2 

ft,2=lih-lh  +  \h  fl/2  =  \h  +  If2  ~  1/3 

l  ft/2  =  ¥ h  -  ffs  +  fh,  {  Zl/2  =  fh  -  5/3  + 1/4. 

Similarly,  for  the  two  points  from  the  right  boundary,  they  are  given  by 


f  fh- 1/2  =  _  f /jv-3  +  t/at-2 

l  fp- 1/2  =  §/iV-3  -  J/jV-2  +  y/tV-1 

[  f V— 1/2  =  \In~2  +  S/jV— 1  +  |/iV, 


{/I  _  13  /  67  x  ,  47  f 

Jn+  1/2  ~  T “  TJ^-3  + 

•4+1/2  =  T/^V-3  -  f /^V-2  + 

4+1/2  =  |/iV-2  -  I/iV-1  +  ^/iV- 


In  all  the  formulas  above,  fj  =  gjUj  are  the  point  values  of  the  flux  function  g(x,Q(x,t))u. 
The  nonlinear  weights  are  defined  by 


wr  = 


ar 


dr 

(e  +  Pr)2 


r  =  1,2,3. 


(2.2) 


Here  dr  are  the  linear  weights  which  can  guarantee  fifth  order  accuracy,  and  £  is  a  small 
positive  number  introduced  to  avoid  the  denominator  to  become  0  and  is  usually  taken  as 
£  =  10-6  in  numerical  tests.  /3r  are  the  so-called  “smoothness  indicators”,  which  measure 
the  smoothness  of  the  function  in  the  relevant  stencils.  We  adopt  the  smoothness  indicators 
introduced  in  [11] 


Pr  = 


1 >-7 . a 

1=1  \aX 


\l  \  2 

Pr(x)  I  dx 


where  pr(x)  is  the  relevant  quadratic  reconstruction  polynomial  to  yield  the  flux  /J+1/2- 
Notice  that,  since  the  coefficient  g(x,Q(x,t ))  is  positive,  upwinding  is  biasing  the  stencil  to 
the  left,  and  the  measurement  of  smoothness  is  for  the  interval  [xj_  1/2,  Xj+1/2]  which  is  to  the 
left  of  the  flux  location  Xj+ 1/2.  We  can  work  out  the  explicit  formulas  for  these  smoothness 
indicators.  These  explicit  formulas,  together  with  the  linear  weights  dr,  are  listed  below.  For 


the  points  inside  the  computational  domain,  with  1  <  j  <  N  —  2,  we  have 

(  Pi  =  if(/l-2  -  fj- 1  +  fj)2  +  iifj-2  -  4/jr-l  +  3 fj)2 

S  P2  =  tK/j-i  -  2 fj  +  fj+i)2  +  \{fj- 1  -  fj+1)2 

{  p3  =  if  {fj  -  2/j+i  +  fj+2)2  +  i(3  fj  -  4/j+i  +  fj+ 2)2 
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and 


di  =  1/10,  d2  =  3/5,  d3  =  3/10. 


For  j  =  N  —  1,  we  have 


Pi  = 


Ip  4-  °±P  4- 

JN- 4  ^  3  JN- 3  ^ 


l  p  _ 

3  J  iV— 4  1  3  •/  jV-3  1  3  JN- 2 

o  _  4  /2  ,25/2  ,10/2 

P2  —  3J N- 3  ■+■  3  J N—2  ~r  3  liV-1 

A  =  l/J-2  +  t/k-1  +  l/iv  -  f /W-2/K-1  +  f/w-2/w 


19  /  /  I  29  /  /  73  /  / 

3"  jN—4jN—3  +  ■3'  J  N—4J  N—2 JJN-3JN-2 

y./'v-3./'v-2  +  y/jv-3/iv-i  —  y/w- 2JW-1 

\fN-2fN  —  y/v-i/v 


and 


di  =  -3/110,  d2  =  47/110,  d3  =  3/5. 


For  j  =  IV,  we  have 


a  _  22  f  2  ,  121  /2  ,  4U  /2  103  /  /  ,  59  /  /  139  /  / 

Pi  —  t/at-4  +  —Jn- 3  +  ~Jn-2 - 3“  J  N —AJ  N —3  +  —  JN-AJN-2  3-JN-3JN-2 

P2  =  ffir- 3  +  ff2N-2  +  ff2N-l-  f./'.V  3«.Y  2  +  f  ,/v  3./.V  ,  -  f /lY-2/lV-l 


ft  =  f/v-2  +  f /v-1  +  f/v  -  VN-2UN-1  +  t/^V-2/^V  ~  ¥/lV-l/lV 


25  /  2 


10  /2 


19 


11 


31 


and 


^  =  3/65,  d2  =  -417/1430  d3  =  137/110. 


Notice  that  some  of  the  linear  weights  for  j  =  IV  —  1  and  j  =  N  are  negative.  We  have 
used  the  technique  introduced  in  [15]  to  treat  these  negative  weights.  We  refer  to  [15]  for 
the  details  and  will  not  describe  them  here. 

For  the  two  points  from  the  left  boundary,  we  could  also  use  the  smoothness  indicators 
and  nonlinear  weights  as  those  for  the  right  boundary,  in  a  mirror  symmetric  fashion.  Our 
computational  experience  however  indicates  that  it  is  better  to  use  simply  the  linear  weights 
for  these  two  points,  which  is  the  practice  that  we  have  adopted  in  the  numerical  experiments 
in  next  section. 


The  global  boundary  condition  at  the  left  is  implemented  by  a  fifth  order  composite  rule 


N 

g0u0  =  C  +  X  '  PjUjAx, 
j=o 


(2.3) 


where  the  special  summation  notation  is  defined  by 


32 

J2'aj 

3=31 


251  299  211  739 

720“*  +  240°*+1  +  24o“*+2  +  72o“*+s  + 


739  211  299 

720°*-3  +  240  “*-2  +  240°*-‘ 


251 

720 


32  4 

+  X  “j 


3=3  1+4 
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if  j‘2  —  ji  >  7.  The  environment  is  computed  also  by  the  same  fifth  order  composite  rule, 
except  for  the  integrals  with  j2  —  ji  <6  and  the  integral  over  the  first  interval  which  is 
computed  avoiding  the  usage  of  u0.  Letting  A  =  (JAuiUi  —  ||cn2w2  +  ^-^3U3  —  ^cn4M4)  Ax  be 
the  fifth  accurate  approximation  to  f*1  u>(x)u(x,  t)  dx,  we  have 


Qj 

Qo 

Q2 

Q3 

Q4 

Qb 

Qe 

Qi 

Qn-i 

Qn-2 


A  +  a  '  UiUiAx  +  '  tUiUiAx,  8  <  j  <  N  —  7, 


A  +  y; '  UjUj Ax,  Qi  =  aA  +  '  UjUjAx,  Qn  =  o:Qo, 

1=1  1=1 

A  (8  5  4  1  \  A, 

aAx  I  -  -a;2u2  +  -u3u3  -  -u4u4  J  +  ^  UiUiAx, 


.  (21  9  15  3  \ 

aAx  I  ycniWi  -  -cu2m2  +  ycn3-u3  ~  g^4«4  J  +  ^i^Ax, 

A  /21  7  29  1  1  \  A, 

«Al  I  y  Witti  -  -^2«2  +  y<^3^3  +  g^4M4  -  —C J5U5  J  +  ^  ^  Ax, 

'  '  i= 4 

A  /21  7  19  17  1  1  \  A, 

«Ax  I  ycniMi  —  -tn2w2  +  —u3u3  +  yCi;4-u4  +  -u>3u3  —  y-u;6-U6  I  +  (UjUjAx, 

'  '  2—5 

A  /21  7  19  2  25  1  1  \ 

aAx  I  ycniMi  —  -U2U2  +  —u3u3  +  -cn4-u4  +  —<^51(5  +  -uJ^ue  —  I 

N 

+  ^  '  uiiUiAx, 


A  /21  7  19  2  25  1  1 

01  Ax  I  y U1U1  —  -cn2w2  +  —u3u3  +  -cn4w4  +  lu3u3  +  —lvqUq  +  -0*77x7  —  — WgMg 

N 

+  '  uiiUiAx, 


/  9  19  5  X 

Ax  (  ~U>nUn  +  —  —UN-2UN-2  +  —  UN-3UN-3 


+a  J  A  +  ^  '  ujiUiAx  j  , 


Ax  (J^UJNUN  +  ^LVn-iUn-!  +  +  a  M  +  '  CUiWjAxj  , 

A  (1  31  7  13  1 

Ax  I  -UjyUpf  +  —Ujst-iUn-i  +  -Un-2UN-2  +  y^V-3MiV-3  —  y UN-4UN-4 


Qn-3 


+a  J  A  +  '  ujiUiAx  j  , 


^  A  (l  31  5  13  1 

Qn-A  —  Ax  -(VnUN  +  —Un-1Un-1  +  ~^N-2UN-2  +  —UN-3UN-3  +  -UN-AUN-A 

V  3  2  4  o  1Z  Z 


:^N-5UN-5 


+  a  I  A  +  ^  '  ujiUiAx  , 


^  .  fl  31  5  25  25 

Qn-5  —  Ax  I  -UisrUjy  +  —  Un-iUn-i  +  ~X>N-2UN-2  +  -^^N-3uN-3  +  ^^N-AUN-A 

1  1  \  /  \ 

+  -<^7V-5WAr-5  ~  —  ^N-SUn-6  J  +  Cl  I  A  +  '  LUiUiAx  |  , 

/ 1  3X  5  25 

Qn-6  =  Ax  I  -UJnUn  +  —UJn-iUn-1  +  -UN-2UN-2  +  —^N-3uN-3  +  ^N-A^N-A 
V  3  24  o  24 


25  1  1  \  /  JV~U  \ 

+  ^^Af-5'uAr-5  +  ^N-6U N-G  —  N-7UN-7  j  +  Cl  I  A  +  '  UiUjAx  I  . 

Notice  that  these  approximations  to  Qj  are  all  hfth  order  accurate.  The  initial  condition  is 


taken  as 


u°j=u0(xj),  j  =  1,2,  •••  ,N. 


Using  the  notation  A  =  £,  we  can  write  the  Euler  forward  time  discretization  of  the  semi¬ 
discrete  scheme  (2.1)  as 

<+1=«”HA(i”1/2-jpi/2)-AtmX.  J  =  1,2. ,JV.  (2.4) 

This,  together  with  the  boundary  condition  (2.3)  for  Uq,  provides  a  complete  description  of 
the  explicit  time  marching  hfth  order  finite  difference  WENO  scheme.  In  order  to  obtain 
higher  order  accuracy  in  time  without  compromising  the  nonlinear  stability  of  the  WENO 
scheme,  we  use  the  total  variation  diminishing  (TVD)  high  order  Runge-Kutta  time  dis¬ 
cretization  in  [18],  see  also  [8,  9].  If  we  denote  the  ordinary  differential  equation  system  (2.1) 


-Uj  —  L(u,  t)j  =  0, 


then  the  third  order  TVD  Runge-Kutta  method  in  [18]  that  we  use  in  next  section  is  given 
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by 


m(1)  =  un  +  AtL(un,tn ) 

w(2)  =  ^n  +  i(«(1)  +  AtL(«(1),tn  +  At))  (2.5) 

■u(2)  +  AfL  ^u(2),t"  + 

Clearly,  this  is  just  three  Euler  forward  time  discretizations  with  suitable  coefficients,  hence 
its  implementation  does  not  introduce  any  difficulty. 


3  Numerical  examples 

In  this  section  we  perform  numerical  experiments  to  demonstrate  the  performance  of  the 
fifth  order  WENO  scheme  developed  in  the  previous  section.  We  use  the  third  order  TVD 
Runge-Kutta  time  discretization  (2.5)  with  a  CFL  condition 

A tn  =  0.6Ax/||^(x,  Qn)  +  m(x,  Qn)  Ax||oo 


unless  otherwise  stated. 

In  our  first  example,  we  use  the  initial  condition  u°{x)  =  —x2  +  x  +  1,  and  take  the 
parameters  and  functions  in  (1.1)  and  (1.2)  as  L  —  1,  a  —  0.5,  u(x)  =  1,  g(x,Q)  = 
(1  —  x)(5  —  x  +  x2 /2  —  Q ),  m{x ,  Q)  =  4  +  2 Q  +  (1  —  x)2/2,  and  / 3{x ,  Q)  —  (1  +  x)(2  —  Q ). 
The  purpose  of  this  example  is  to  show  that  the  WENO  scheme  is  non-oscillatory  in  the 
presence  of  solution  discontinuities.  For  this  purpose  we  take  C{t)  =  3,  which  causes  an 
incompatibility  of  the  boundary  data  and  the  initial  condition  at  the  origin  and  generates 
a  discontinuity  in  the  solution  which  moves  from  the  left  boundary  into  the  computational 
domain.  The  numerical  solutions  using  N  =  100  uniformly  spaced  grid  points  for  the  first 
order  monotone  scheme  [14],  the  second  order  high  resolution  scheme  [14],  and  the  fifth  order 
WENO  scheme  are  plotted  in  Fig.  3.1.  We  can  see  clearly  that  our  fifth  order  WENO  scheme 
can  resolve  the  discontinuity  sharper  than  the  Erst  order  monotone  scheme  and  the  second 
order  high  resolution  scheme,  without  introducing  spurious  numerical  oscillations. 
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Figure  3.1:  Numerical  solutions  with  N  =  100  uniform  grid  points,  using  the  first  order 
monotone  scheme  in  [14]  (triangles),  the  second  order  high  resolution  scheme  in  [14]  (circles), 
and  the  fifth  order  WENO  scheme  (solid  line). 

The  second  example  is  chosen  to  demonstrate  that  our  fifth  order  WENO  scheme  can 
achieve  its  designed  accuracy  for  smooth  solutions.  For  this  purpose  we  take  g(x,  Q )  = 
— 2ex~1  +  2,  / 3(x,Q )  =  2,  m(x,Q)  —  1,  L  —  1,  u(x)  —  1,  a  —  0.5,  the  initial  condition 
u°(x)  =  e~x,  and  C(t)  =  0.  For  these  choices  the  equation  (1.1)  has  a  smooth  exact  solution 
u(x,t)  =  et~x .  The  L°°  and  L 1  errors  and  orders  of  accuracy  of  our  fifth  order  WENO 
scheme  are  listed  in  Table  3.1.  We  can  see  that  the  designed  order  of  accuracy  is  obtained 
for  this  smooth  solution  in  the  L 1  norm.  The  order  in  the  L°°  norm  is  reduced  by  one  order 
because  near  the  boundary,  the  cancellation  of  the  neighboring  numerical  fluxes  to  generate 
an  extra  order  of  accuracy  is  no  longer  valid.  For  this  test  case  we  further  reduce  the  time 
step  to  ensure  that  the  spatial  error  dominates  in  the  solution. 

Our  third  example  is  from  computational  biology  [2] ,  for  the  simulation  of  the  population 
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Table  3.1:  L°°  and  L 1  errors  and  numerical  orders  of  accuracy  of  the  fifth  order  WENO 
scheme  using  N  uniformly  spaced  mesh  points. 


N 

L°°  error 

order 

L 1  error 

order 

10 

0.21E-03 

0.25E-04 

20 

0.95E-05 

4.46 

0.64E-06 

5.28 

40 

0.52E-06 

4.20 

0.19E-07 

5.07 

80 

0.30E-07 

4.10 

0.59E-09 

5.03 

160 

0.18E-08 

4.08 

0.18E-10 

5.02 

320 

0.10E-09 

4.12 

0.49E-12 

5.20 

evolution  of  Gambussia  affinis.  This  model  is  given  by  the  equation  (1.1),  written  in  a  slightly 
different  form  as 


ut  +  (g(x,  t )  u)x  +  m(x,  Q(t),  t)u  —  0,  (x,  t )  G  [9,  63]  x  (0,  T } 


-63 


g(9,t)u(9,t)  —  /  P(x,t)  u(x,t)dx,  f  G  (0,  T] 
J  9 

u(x,  0)  =  u°(x),  x  G  [9,63] 


(3.1) 


The  non-local  term  Q(t)  is  defined  by 


r-63 


Q(t)  —  /  c j(x)u(x,t)dx 


(3.2) 


The  initial  condition  and  the  functions  in  (3.1)  and  (3.2)  are  defined  as  in  [2],  In  particular, 

P(x,  t)  =  P(x)Tp(t),  g(x,  t)  =  g(x)Tg(t),  m(x,  Q(t),t )  =  m(x,  Q(t))Tm(t), 

where  /3(x)  is  a  smooth  spline  function  to  fit  the  data  in  Krumholtz  [13]  by  using  a  MATLAB 
function  called  esaps  (see  Fig.  1  in  [2]),  Tp(t )  is  defined  by 


_  t- 30  _|_  h-30)2  ■ 


Mt)  = 


o. 


't-120\3/'-i  I  t— 90  ,  (t— 90)2  \ 
30  /  V1  “r  10  "l“  150  b 


0  <  t  <  30 
30  <  t  <  90 
90  <  t  <  120 
120  <  t  <  365 


and  it  is  periodically  extended  thereafter 


Tp(t  +  365 n)  =  Tp(t),  n  =  1,  2, 
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The  function  g(x)  is  defined  as 


g(x) 


9  <  x  <  63 


and  the  function  Tg(t)  is  defined  as 


Tg(t)  =  0.2 +  0.8  Tp{t). 


The  function  m(x,  Q )  is  given  by 


m(x,  Q ) 


0.1  exp(— C/Q),  9  <  x  <  31 

0.1  exp  (—C/Q)  —  (0.023  —  0.1  exp(— C/Q)) 

x(x-  31)3(1  -  3(x  -  32) (65  -  2x)),  31  <  x  <  32 

0.023,  32  <  x  <  63 


where  the  constant  C  will  be  prescribed  later.  The  function  Tm(t )  is  given  by 


(3.3) 


Tm{t)  =  2  —  T@(t), 


and  finally  the  function  u(x)  is  given  by 


oj(x) 


2,  9  <  x  <  30 

-2(x  -  31)3(1  +  3(x  -  30) (2a:  -  59)),  30  <  x  <  31 

0.  31  <  x  <  63 


The  initial  condition  is  given  as 


Uq(x)  =  { 


0, 


5(1  +  ^)3, 

5  +  15(^)  +  15(^)2  +  30( 
5(2  -^)3, 

0. 


x—38 

4 


x— 46  \ 
4 


9  <  x  <  34 
34  <  x  <  38 
38  <  x  <  42 
42  <  x  <  46 
46  <  x  <  63 


We  note  that  not  all  the  assumptions  (H1)-(H6)  outlined  in  the  introduction  are  satisfied  in 
this  example.  However,  our  fifth  order  WENO  scheme  performs  nicely  and  gives  accurate 
results  with  far  fewer  grid  points  than  the  second  order  schemes  in  [2]  and  [14].  In  Fig.  3.2, 
we  plot  the  population  density  u  at  t  —  365  using  the  fifth  order  WENO  scheme  with 


N  =  135  uniformly  spaced  grid  points.  In  this  simulation  the  constant  C  in  (3.3)  is  taken  as 
2000  as  in  [2],  For  the  purpose  of  comparison,  the  simulation  result  using  the  second  order 
high  resolution  scheme  in  [14]  with  N  =  540  uniformly  spaced  grid  points  is  also  plotted  in 
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Fig.  3.2.  We  can  observe  that  the  WENO  scheme  gives  better  resolution  than  the  second 
order  scheme  even  though  the  latter  uses  fourfold  as  many  grid  points.  Comparing  with  the 
second  order  box  method  using  N  =  730  grid  points  in  [2],  Fig.  4,  we  can  see  that  the  second 
order  box  method  in  [2]  performs  similarly  as  the  second  order  high  resolution  scheme  in 
[14],  and  the  fifth  order  WENO  scheme  clearly  outperforms  both  of  them. 


Figure  3.2:  Population  density  u  as  a  function  of  the  length  at  t  —  365.  Solid  line:  solution 
of  the  fifth  order  WENO  scheme  with  N  =  135  uniformly  spaced  grid  points;  Circle  symbols: 
solution  of  the  second  order  high  resolution  scheme  [14]  with  N  =  540  uniformly  spaced  grid 
points. 


In  order  to  fully  assess  the  resolution  power  of  the  high  order  WENO  scheme  for  long 
time  simulation  with  coarse  meshes,  we  simulate  this  model  for  10  years  and  plot  the  total 
population  (the  integral  of  the  density  u  over  the  length)  in  Fig.  3.3,  for  two  different  values 
of  C  in  (3.3),  namely  C  =  2000  for  the  left  picture  and  C  =  200000  for  the  right  picture, 
as  in  [2],  Figures  7  and  12.  For  this  test  problem,  we  perform  numerical  tests  with  different 
number  of  mesh  points  and  report  the  results  using  the  coarsest  meshes  of  different  numerical 
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methods  to  obtain  the  visually  satisfactory  resolution.  We  are  pleasantly  surprised  to  observe 
that  the  WENO  scheme  with  only  N  —  20  uniformly  spaced  grid  points  is  enough  to  yield 
the  satisfactory  resolution  as  shown  in  Fig.  3.3,  while  the  second  order  high  resolution  scheme 
in  [14]  needs  N  =  108  points  to  achieve  comparable  resolution. 


Figure  3.3:  Evolution  of  the  total  population  for  ten  years.  Left:  C  =  2000;  Right:  C  = 
200000.  Solid  line:  solution  of  the  fifth  order  WENO  scheme  with  N  =  20  uniformly  spaced 
grid  points;  Circle  symbols:  solution  of  the  second  order  high  resolution  scheme  [14]  with 
N  =  108  uniformly  spaced  grid  points. 


4  Concluding  remarks 


We  have  developed  a  fifth  order  explicit  finite  difference  WENO  scheme  for  solving  a  hierar¬ 
chical  size-structured  population  model  with  nonlinear  growth,  mortality  and  reproduction 
rates,  which  contains  global  terms  both  for  the  boundary  condition  and  for  the  coefficients  in 
the  equations.  Numerical  results  are  provided  to  demonstrate  the  capability  of  this  WENO 
scheme  in  resolving  smooth  as  well  as  discontinuous  solutions.  For  smooth  solutions  the  fifth 
order  WENO  scheme  achieves  its  designed  order  of  accuracy.  For  discontinuous  solutions, 
the  WENO  scheme  gives  sharp  and  non-oscillatory  discontinuity  transitions.  An  application 
of  the  scheme  to  a  problem  in  computational  biology  for  the  evolution  of  the  population  of 
Gambussia  affinis  indicates  that  the  WENO  scheme  can  achieve  good  resolution  for  long 
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time  simulation  with  very  coarse  meshes,  hence  it  is  much  more  efficient  than  lower  order 
schemes  for  such  simulations. 
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